package org.farrell.EOSCalc;

import static java.lang.Math.abs;
import static java.lang.Math.log;
import static java.lang.Math.sqrt;

public class PREOS extends CubicEOS implements Compound {
	// This implements the Peng-Robinson equation of state

	public static double kappa;
	
	public PREOS(double criticalTemperature, double criticalPressure, double acentricFactor) {
		super(criticalTemperature, criticalPressure, acentricFactor);
		kappa = 0.37464+w*(1.54226-w*0.26993);
	}
	
	public double getAlpha(double Temperature){
		double Tr = Temperature/Tc;
		double inter = 1+kappa*(1-sqrt(Tr));
		return inter*inter;
	}
	
	@Override
	public double getA(double Temperature, double Pressure){
		// Calculates the A constant for the Peng-Robinson EOS
		double alpha = getAlpha(Temperature);
		double invTr = Tc/Temperature;
		double Pr = Pressure/Pc;
        return 0.45723553*alpha*Pr*invTr*invTr;
	}
	
	@Override
	public double getB(double Temperature, double Pressure){
		// Calculates the B constant for the Peng-Robinson EOS
		double invTr = Tc/Temperature;
		double Pr = Pressure/Pc;
        return 0.07779607*Pr*invTr;
	}

	@Override
	public Double[] getZ(double a, double b){
		// Calculates the roots of the EOS
		// returns a double of size numRoots		
	    double a2 = b-1;
	    double a1 = a - b*(3*b+2);
	    double a0 = b*(b*(1+b)-a);
	    return cubicRoots(a0,a1,a2);
	}
	
	@Override
	public double getHdepart(double z, double a, double b, double temp){
		double m = Tc*getAlpha(temp)/temp;
    	return z-1-2.0779606467423424*(m+kappa*sqrt(m))*log(abs((z+2.4142135623730949*b)/(z-0.4142135623730949*b)));
	}
	
	@Override
	public double getSdepart(double z, double a, double b, double temp){
		double m = Tc*getAlpha(temp)/temp;
    	return log(abs(z-b))-2.0779606467423424*kappa*sqrt(m)*log(abs((z+2.4142135623730949*b)/(z-0.4142135623730949*b)));
	}
	
	@Override
	public double getGdepart(double z, double a, double b, double temp){
		return z-1-log(abs(z-b))-0.35355339059327373*(a/b)*log(abs((z+2.4142135623730949*b)/(z-0.4142135623730949*b)));
	}
    

}
